inline
dct2d::dct2d(const  uword  in_N)
: N(in_N)
{
}


inline 
mat
dct2d::get_coeff_matrix(const mat in_frame)
{
  frame = in_frame;
  uword ro =  frame.n_rows;
  uword co =  frame.n_cols;
  mean_frame  = mean(mean(frame));
  frame.reshape(ro*co,1);
  std_frame=as_scalar(stddev(frame));
  frame.reshape(ro,co);
  
  uword k = 0;
  Cuv=zeros(ro,co);
  
  double lw = (co - N) + 1;
  double lh = (ro - N) + 1;
  double l  = lw*lh; 
  
  
  mat coef_matrix = zeros(N*N, l);
  //mat coef_matrizTopLeft = zeros(15,64); 
  mat coef;
  //mat coef2;
  //vec tmp;
  
      
  for(uword row = 0; row < ro - N; row+=N)
  {
    uword r_in = row;
    uword r_fi = row + (N-1);
    
    for(uword col = 0; col < co - N; col+=N)
    {
      uword c_in = col;
      uword c_fi = col + (N-1);
           
       
      coef  = calc_dct_one_block ( r_in, r_fi, c_in, c_fi );
      
      Cuv.submat(r_in,c_in,r_fi,c_fi) = coef; //Using all features
      
      coef.reshape(N*N,1);
      coef_matrix.col(k)=coef;
      k++;
    }
  }
  
  return coef_matrix;

}

inline
mat 
dct2d::calc_dct_one_block (const uword ri, const uword rf, const uword ci, const uword cf)
{
  mat subframe = frame.submat(ri,ci,rf,cf);
  //int x, y, u, v;
  mat dct = zeros(N,N);
  double s;
  
  for (uword u = 0; u < N; u++) 
    for (uword v = 0; v < N; v++)
    {
      s = 0;
      for (uword x = 0; x < N; x++)
	for (uword y = 0; y < N; y++)
	{
	  double Cu = ((u == 0) ? 1 / sqrt(2) : 1);
	  double Cv = ((v == 0) ? 1 / sqrt(2) : 1);
	  //double subframe_xy = (subframe(x,y) - mean_frame)/std_frame; 
	  double subframe_xy = subframe(x,y);
	  s += subframe_xy * cos((2 * x + 1) * u * datum::pi / (2*N)) *
	  cos((2 * y + 1) * v * datum::pi / (2*N)) * Cu * Cv ;
	}
	dct(u,v) = s*2/N; // = s*2/N
	//cout << dct(u,v) << "..";
    }
    return dct;
}

inline
mat 
dct2d::calc_dct_one_block16 (const uword ri, const uword rf, const uword ci, const uword cf) //
{
  mat subframe = frame.submat(ri,ci,rf,cf);
  //int x, y, u, v;
  mat dct = zeros(4,4);
  double s;
  
  for (uword u = 0; u < 4; u++) // using the first 16 function bases (0-3)*(0-3)
    for (uword v = 0; v < 4; v++) // using the first 16 function bases
    {
      s = 0;
      for (uword x = 0; x < N; x++)
	for (uword y = 0; y < N; y++)
	{
	  double Cu = ((u == 0) ? 1 / sqrt(2) : 1);
	  double Cv = ((v == 0) ? 1 / sqrt(2) : 1);
	  //double subframe_xy = (subframe(x,y) -mean_frame)/std_frame; 
	  double subframe_xy = subframe(x,y);
	  s += subframe_xy * cos((2 * x + 1) * u * datum::pi / (2*N)) *
	  cos((2 * y + 1) * v * datum::pi / (2*N)) * Cu * Cv ;
	}
	dct(u,v) = s*2/N; // = s*2/N
    }
    return dct;
}




inline 
mat
dct2d::get_coeff_overlap_7(const mat in_frame)//Using whole feature set, overlapping 7 pixels, jumping 1
{
  frame = in_frame;
  uword ro =  frame.n_rows;
  uword co =  frame.n_cols;
  
  mean_frame  = mean(mean(frame));
  frame.reshape(ro*co,1);
  std_frame=as_scalar(stddev(frame));
  frame.reshape(ro,co);
  
  uword k = 0;
  mat coef;
  
  double lw = (co - N) + 1;
  double lh = (ro - N) + 1;
  double l  = lw*lh; 
   
   
  //uword l = pow( (N*N-N+1), 2);
  
  
  mat coef_matrix = zeros(N*N,l);
  
  for(uword row = 0; row <= ro-N; ++row)
  {
    uword r_in = row;
    uword r_fi = row + (N-1);
    
    for(uword col = 0; col <= co-N; ++col)
    {
      uword c_in = col;
      uword c_fi = col + (N-1);
      coef  = calc_dct_one_block( r_in, r_fi, c_in, c_fi );
      coef.reshape(N*N,1);
      coef_matrix.col(k) = coef;
      k++;
    }
  }
  return coef_matrix;
}

inline 
mat
dct2d::get_coeff_overlap_6(const mat in_frame)//Using Dim = 15 feature set, overlapping 6 pixels JUMPING 2
{
  frame 	= in_frame;
  uword ro	= frame.n_rows;
  uword co	= frame.n_cols;
  
  mean_frame  = mean(mean(frame));
  frame.reshape(ro*co,1);
  
  std_frame   = as_scalar(stddev(frame));
  frame.reshape(ro,co);
  
  mat coef;
  uword k    = 0;
  
  uword jump = 2;
  
  
  double lw = floor((co - N)/jump + 1);
  double lh = floor((ro - N)/jump + 1);
  double l  = lw*lh; 
   //cout << "l= "<< l << endl;
   //cout << "Actual size = ( " << frame.n_rows << " , " << frame.n_cols  << " )"<< endl;

  //uword l    = pow( (w-N)/jump + 1, 2); //if h = w (usually 64*64) 
  
  mat coef_matrix = zeros (15,l); //using features 1-15
  //uword last_nw = w-N;
  
  
  for(uword row = 0; row <= ro - N ; row = row + jump)
  {
    uword r_in = row;
    uword r_fi = row + (N-1);
    
    for(uword col = 0; col <= co - N ; col = col + jump)
    {
      uword c_in = col;
      uword c_fi = col + (N-1);
      coef  = calc_dct_one_block16( r_in, r_fi, c_in, c_fi );
      //cout << "coef.n_elem= " <<coef.n_elem << endl;
      coef.reshape(16,1);
      vec tmp = coef;
      //cout << "features per block:" << endl  << tmp.t() << endl;
      coef_matrix.col(k) = tmp.subvec(1,15);
     k++;
    }
  }
  //cout << k << endl;
  //cout << l << endl;
  return coef_matrix;
}



// reconst_frame is only useful when all features are calculated (get_coeff_matrix)
inline 
mat
dct2d::reconst_frame()
{
  mat new_frame=zeros(frame.n_rows, frame.n_cols );
  mat tmp;
  
  for(uword row = 0; row < frame.n_rows - N; row+=N)
  {
    uword r_in = row;
    uword r_fi = row + (N-1);
    
    for(uword col = 0; col < frame.n_cols -N; col+=N)
    {
      uword c_in = col;
      uword c_fi = col + (N-1);
      tmp  = calc_Idct_one_block (r_in, r_fi, c_in, c_fi );
      new_frame.submat(r_in,c_in,r_fi,c_fi) = tmp; //Using all features
    }
  }
  return new_frame;
}


inline
mat
dct2d::calc_Idct_one_block (const uword ri, const uword rf, const uword ci, const uword cf)
{
  int i, j, u, v;
  double s;
  mat out = zeros(N,N);
  mat dct = Cuv.submat(ri,ci,rf,cf);
  
  for (i = 0; i < N; i++)
    for (j = 0; j < N; j++)
    {
      s = 0;
      for (u = 0; u < N; u++)
	for (v = 0; v < N; v++)
	{
	  double Cu = ((u == 0) ? 1 / sqrt(2) : 1);
	  double Cv = ((v == 0) ? 1 / sqrt(2) : 1);
	  
	  s += dct(u,v) * cos((2 * i + 1) * u * datum::pi/ (2*N)) *
	  cos((2 * j + 1) * v * M_PI / (2*N)) * Cu * Cv;
	}
	out(i,j) = s*2/N;
    }
    //cout << "Idct goodiz " << endl;
    return out;
}

